module moments3
!moments of trivariate Mardia's normality test
implicit none
contains
!!!!
subroutine mardiamvu(MD,VD,matA,matB,b,eta,si,sigma,dmut,dvsig,pth)
use linear_operators
use testmoments
use testmoments2
use testmomsrot
use matrixfuns
integer, parameter:: n=3,p=10
integer, intent(in)::pth
double precision, intent(in)::b(n),eta,si,dmut(pth,n),dvsig(pth,n*n),sigma(n,n)
double precision, intent(out)::matA(pth,pth),matB(p+1,pth),MD(pth+p+1),VD(pth+p+1,pth+p+1)
double precision:: isigma(n,n),rsigma(n,n),irsigma(n,n),beta(n),dmach,tol&
&,matthmsk(n,p),matthvsk(n*n,p),msk(p,1),m32,m14,m122&
&,matV(n,n),matL(n*(n+1)/2,n*n),matD(n*n,n*(n+1)/2),matDm(n*(n+1)/2,n*n),mstar1(n,p),mstar2(n*n,p)&
&,dep(n,pth),dep2(n,pth),me(n,1),m3,m12,irsigmah(n,n),drvsig(n*n,n*n)
integer i,j,irank

isigma=.i.sigma
tol = 100.0*dmach(4)
call dchfac(n,sigma,n,tol,irank,rsigma,n)
rsigma=.t.rsigma
irsigma=.i.rsigma
beta=matmul(.t.rsigma,b)

if (dot_product(beta(2:n),beta(2:n))<1.0d-10) then
	matV=eye(n)
else
!	Sigma=eye(n)
   matV(1,:)=(/-0.57735026918963d0, 0.81649658092773d0,0.0d0/)
   matV(2,:)=(/-0.57735026918963d0,-0.40824829046386d0,-0.70710678118655d0/)
   matV(3,:)=(/-0.57735026918963d0,-0.40824829046386d0, 0.70710678118655d0/)

!	matV(:,1)=beta/dsqrt(dot_product(beta,beta))
!	matV(1,2:3)=(/-0.50497514502648d0,-0.58533678384716d0/)
!	matV(2,2:3)=(/0.78236638609871d0,-0.0d0/)
!	matV(3,2:3)=(/-0.36455855607615d0,0.81079026232156d0/)

!	print*, "Please, supply orthogonal matrix!"
!	print*, "beta", beta
!	stop
	!matV(:,1)=beta/dsqrt(dot_product(beta,beta))
end if
irsigmah=matV.tx.irsigma

beta(1)=dsqrt(dot_product(beta,beta))
beta(2:3)=(/0.0d0,0.0d0/)

forall(i=1:pth)
MD(i)=0.0d0
end forall
call meanku(MD(pth+1),beta,eta,si)
call meansk(MD(pth+2:pth+p+1),beta,eta,si)

VD(1:pth,1:pth)=matmul(dmut,matmul(isigma,.t.dmut))&
&+.25d0*(((dvsig.x.kron(.t.irsigmah,.t.irsigmah))&
&.x.(emempkememp(beta,eta,si)-(vec(dble(eye(n))).xt.vec(dble(eye(n))))))&
&.x.(kron(irsigmah,irsigmah).xt.dvsig))&
&+.5d0*((dvsig.xt.kron(irsigmah,irsigmah)).x.(emkememp(beta,eta,si).x.(irsigmah.xt.dmut)))&
&+.5d0*(.t.(((dvsig.xt.kron(irsigmah,irsigmah)).x.(emkememp(beta,eta,si).x.(irsigmah.xt.dmut)))))

call varku(VD(pth+1,pth+1),beta,eta,si)
call varsk(VD(pth+2:pth+p+1,pth+2:pth+p+1),beta,eta,si)
call covkusk(VD(pth+1,pth+2:pth+p+1),beta,eta,si)

VD(1:pth,pth+1)=matmul(dmut.xt.irsigmah,ejm2em(beta,eta,si))&
&+.5d0*matmul(dvsig.xt.kron(irsigmah,irsigmah),vec2(ejm2ememp(beta,eta,si)-ejm2(beta,eta,si)*eye(n)))

forall(i=1:n,j=1:p)
	matthmsk(i,j)=0.0d0
end forall
matthmsk(1,1)=eubeta4(dsqrt(dot_product(beta,beta)),eta,si)
matthmsk(2,2)=ehr(-1.0d0/(2.0d0*eta),(1.0d0-si)/si,2)*3.0d0
matthmsk(3,3)=matthmsk(2,2)
matthmsk(2,4)=eub2hn(dsqrt(dot_product(beta,beta)),eta,si,1)
matthmsk(3,5)=matthmsk(2,4)
matthmsk(1,6)=matthmsk(2,4)
matthmsk(3,7)=ehr(-1.0d0/(2.0d0*eta),(1.0d0-si)/si,2)
matthmsk(1,8)=matthmsk(2,4)
matthmsk(2,9)=ehr(-1.0d0/(2.0d0*eta),(1.0d0-si)/si,2)
forall(i=1:n*n,j=1:p)
	matthvsk(i,j)=0.0d0
end forall
m32=eub3hn(dsqrt(dot_product(beta,beta)),eta,si,1)
m14=eub1hn(dsqrt(dot_product(beta,beta)),eta,si,2)*3.0d0
m122=eub1hn(dsqrt(dot_product(beta,beta)),eta,si,2)

matthvsk(1,1)=eubeta5(dsqrt(dot_product(beta,beta)),eta,si)
matthvsk(1,6)=m32
matthvsk(1,8)=m32

matthvsk(2,2)=m14
matthvsk(2,4)=m32
matthvsk(2,9)=m122

matthvsk(3,3)=m14
matthvsk(3,5)=m32
matthvsk(3,7)=m122

matthvsk(4,2)=m14
matthvsk(4,4)=m32
matthvsk(4,9)=m122

matthvsk(5,1)=m32
matthvsk(5,6)=m14
matthvsk(5,8)=m122

matthvsk(6,10)=m122

matthvsk(7,3)=m14
matthvsk(7,5)=m32
matthvsk(7,7)=m122

matthvsk(8,10)=m122

matthvsk(9,1)=m32
matthvsk(9,6)=m122
matthvsk(9,8)=m14
msk(:,1)=MD(pth+2:pth+p+1)



VD(1:pth,pth+2:pth+p+1)=matmul(dmut.xt.irsigmah,matthmsk)&
&+.5d0*matmul(dvsig.xt.kron(irsigmah,irsigmah),matthvsk-(vec(dble(eye(n))).xt.msk))

forall(i=1:pth,j=pth+1:pth+p+1)
VD(j,i)=VD(i,j)
end forall
VD(pth+2:pth+p+1,pth+1)=VD(pth+1,pth+2:pth+p+1)


matA=matmul(dmut,matmul(isigma,.t.dmut))&
&+.5d0*matmul(dvsig,matmul(kron(isigma,isigma),.t.dvsig))

matB(1,:)=-4.0d0*matmul(dmut.xt.irsigmah,ejmem(beta,eta,si))&
&-2.0d0*matmul(dvsig.xt.kron(irsigmah,irsigmah),vec2(ejmememp(beta,eta,si)))

forall(i=1:n*(n+1)/2,j=1:n*n)
	matL(i,j)=0.0d0
	matD(j,i)=0.0d0
end forall
matL(1:3,1:3)=eye(3)
matL(4:5,5:6)=eye(2)
matL(6,9)=1.0d0

matD(1:3,1:3)=eye(3)
matD(4,2)=1.0d0
matD(7,3)=1.0d0
matD(5:6,4:5)=eye(2)
matD(8,5)=1.0d0
matD(9,6)=1.0d0
matDm=matmul(.i.(matD.tx.matD),.t.matD)
mstar1=-(matmul(.t.matV,irsigma).xt.dmut)

drvsig=.5d0*((matL.tx.(.i.((matDm.x.kron(rsigma,dble(eye(n)))).xt.matL))).x.matDm)
mstar2=-kron(.t.matV,(.t.matV).x.irsigma)&
&.x.(drvsig.xt.dvsig)



m3=eubeta3(dsqrt(dot_product(beta,beta)),eta,si)
m12=eub1hn(dsqrt(dot_product(beta,beta)),eta,si,1)
!1
me(:,1)=(/m3,0.0d0,0.0d0/)
dep=(mstar1+(kron(me,dble(eye(n))).tx.mstar2))
matB(2,:)=3.0d0*dep(1,:)
!2
me(:,1)=(/m12,0.0d0,0.0d0/)
dep=(mstar1+(kron(me,dble(eye(n))).tx.mstar2))
matB(3,:)=3.0d0*dep(2,:)
!3
me(:,1)=(/m12,0.0d0,0.0d0/)
dep=(mstar1+(kron(me,dble(eye(n))).tx.mstar2))
matB(4,:)=3.0d0*dep(3,:)
!4
me(:,1)=(/0.0d0,m12,0.0d0/)
dep=(kron(me,dble(eye(n))).x.mstar2)
me(:,1)=(/m3,0.0d0,0.0d0/)
dep2=(mstar1+(kron(me,dble(eye(n))).tx.mstar2))
matB(5,:)=2.0d0*dep(1,:)+dep2(2,:)
!5
me(:,1)=(/0.0d0,0.0d0,m12/)
dep=(kron(me,dble(eye(n))).x.mstar2)
me(:,1)=(/m3,0.0d0,0.0d0/)
dep2=(mstar1+(kron(me,dble(eye(n))).tx.mstar2))
matB(6,:)=2.0d0*dep(1,:)+dep2(3,:)
!6
me(:,1)=(/0.0d0,m12,0.0d0/)
dep=(kron(me,dble(eye(n))).x.mstar2)
me(:,1)=(/m12,0.0d0,0.0d0/)
dep2=(mstar1+(kron(me,dble(eye(n))).tx.mstar2))
matB(7,:)=2.0d0*dep(2,:)+dep2(1,:)
!7
me(:,1)=(/m12,0.0d0,0.0d0/)
dep2=(mstar1+(kron(me,dble(eye(n))).tx.mstar2))
matB(8,:)=dep2(3,:)
!8
me(:,1)=(/0.0d0,0.0d0,m12/)
dep=(kron(me,dble(eye(n))).x.mstar2)
me(:,1)=(/m12,0.0d0,0.0d0/)
dep2=(mstar1+(kron(me,dble(eye(n))).tx.mstar2))
matB(9,:)=2.0d0*dep(3,:)+dep2(1,:)
!9
me(:,1)=(/m12,0.0d0,0.0d0/)
dep2=(mstar1+(kron(me,dble(eye(n))).tx.mstar2))
matB(10,:)=dep2(2,:)
!10
me(:,1)=(/0.0d0,0.0d0,m12/)
dep=(kron(me,dble(eye(n))).tx.mstar2)
me(:,1)=(/0.0d0,m12,0.0d0/)
dep2=((kron(me,dble(eye(n))).tx.mstar2))
matB(11,:)=dep(2,:)+dep2(3,:)


end subroutine mardiamvu

subroutine mardiamv(MD,VD,b,eta,si,sigma)
use linear_operators
integer, parameter:: n=3,p=10
double precision, intent(in)::b(n),eta,si,sigma(n,n)
double precision, intent(out)::MD(p+1),VD(p+1,p+1)
double precision:: isigma(n,n),rsigma(n,n),irsigma(n,n),beta(n),dmach,tol
integer i,j,irank

isigma=.i.sigma
tol = 100.0*dmach(4)
call dchfac(n,sigma,n,tol,irank,rsigma,n)
rsigma=.t.rsigma
irsigma=.i.rsigma
beta=matmul(.t.rsigma,b)

call meanku(MD(1),beta,eta,si)
call meansk(MD(2:p+1),beta,eta,si)

call varku(VD(1,1),beta,eta,si)
call varsk(VD(2:p+1,2:p+1),beta,eta,si)
call covkusk(VD(1,2:p+1),beta,eta,si)
VD(2:p+1,1)=VD(1,2:p+1)
end subroutine mardiamv
!!!!
subroutine meansk(matmsk,beta,eta,si)
!Mean of the moment conditions of the skewness comp.
use ghs
use testmomsrot
double precision matmsk(10),beta(3),eta,si,beta1
integer i

forall(i=1:10)
	matmsk(i)=0.0d0
end forall
beta1=dsqrt(dot_product(beta,beta))

matmsk(1)=eubeta3(beta1,eta,si)
matmsk(6)=eub1hn(beta1,eta,si,1)
matmsk(8)=eub1hn(beta1,eta,si,1)
end subroutine meansk

subroutine meanku(matmku,beta,eta,si)
!Mean of the moment conditions of the kurtosis comp.
use ghs
use testmoments
double precision matmku,beta(3),betarot(3),eta,si
integer i
betarot(1)=dsqrt(dot_product(beta,beta))
betarot(2:3)=(/0.0d0,0.0d0/)
matmku=ejm2(betarot,eta,si)-3.0d0*5.0d0
end subroutine meanku
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine varsk(matvsk,beta,eta,si)
!Mean of the moment conditions of the skewness comp.
use ghs
use testmomsrot
use linear_operators
double precision matvsk(10,10),beta(3),eta,si,beta1,matmsk(10,1),m42&
&,m24,nu,gam,m222
integer i,j

forall(i=1:10,j=1:10)
	matvsk(i,j)=0.0d0
end forall
beta1=dsqrt(dot_product(beta,beta))
nu=-1.0d0/(2.0d0*eta)
gam=(1.0d0-si)/si
call meansk(matmsk(:,1),beta,eta,si)

m42=eub4hn(beta1,eta,si,1)			!E(e1^4*ej^2)
m222=eub2hn(beta1,eta,si,2)			!E(e1^2*ej^2*ek^2)
m24=eub2hn(beta1,eta,si,2)*3.0d0    !E(e1^2*ej^4)
matvsk(1,1)=eubeta6(beta1,eta,si)
matvsk(6,1)=m42
matvsk(1,6)=m42
matvsk(8,1)=m42
matvsk(1,8)=m42

matvsk(2,2)=ehr(nu,gam,3)*15.0d0
matvsk(4,2)=m24
matvsk(2,4)=m24
matvsk(9,2)=ehr(nu,gam,3)*3.0d0
matvsk(2,9)=matvsk(9,2)

matvsk(3,3)=matvsk(2,2)
matvsk(5,3)=m24
matvsk(3,5)=m24
matvsk(7,3)=ehr(nu,gam,3)*3.0d0
matvsk(3,7)=matvsk(7,3)

matvsk(4,4)=m42
matvsk(9,4)=m222
matvsk(4,9)=m222

matvsk(5,5)=m42
matvsk(7,5)=m222
matvsk(5,7)=m222

matvsk(6,6)=m24
matvsk(8,6)=m222
matvsk(6,8)=m222

matvsk(7,7)=ehr(nu,gam,3)*3.0d0
matvsk(8,8)=m24
matvsk(9,9)=ehr(nu,gam,3)*3.0d0
matvsk(10,10)=m222


matvsk=matvsk-(matmsk.xt.matmsk)
end subroutine varsk

subroutine varku(matvku,beta,eta,si)
!Mean of the moment conditions of the kurtosis comp.
use ghs
use testmoments
double precision, parameter:: n=3.0d0
double precision matvku,beta(3),betarot(3),eta,si,matmku
integer i
betarot(1)=dsqrt(dot_product(beta,beta))
betarot(2:3)=(/0.0d0,0.0d0/)
matmku=ejm2(betarot,eta,si)-3.0d0*5.0d0
matvku=ejm4(betarot,eta,si)+n*n*(n+2.0d0)*(n+2.0d0)&
&-2.0d0*n*(n+2.0d0)*ejm2(betarot,eta,si)-matmku*matmku
end subroutine varku
!!!

subroutine covkusk(matcsk,beta,eta,si)
!Mean of the moment conditions of the skewness comp.
use ghs
use testmomsrot
double precision matcsk(10),matmsk(10),matmku,beta(3),eta,si,beta1
integer i

beta1=dsqrt(dot_product(beta,beta))
call meansk(matmsk,beta,eta,si)
call meanku(matmku,beta,eta,si)
forall(i=1:10)
matcsk(i)=0.0d0
end forall
matcsk(1)=eubeta7(beta1,eta,si)+4.0d0*eub5hn(beta1,eta,si,1)+8.0d0*eub3hn(beta1,eta,si,2)
matcsk(6)=eub5hn(beta1,eta,si,1)+8.0d0*eub3hn(beta1,eta,si,2)+24.0d0*eub1hn(beta1,eta,si,3)
matcsk(8)=eub5hn(beta1,eta,si,1)+8.0d0*eub3hn(beta1,eta,si,2)+24.0d0*eub1hn(beta1,eta,si,3)
matcsk=matcsk-3.0d0*5.0d0*matmsk-matmku*matmsk
end subroutine covkusk

!!!
subroutine var0(wmat)
!Variance of the moment conditions under the null
use ghs
use testmomsrot
use linear_operators
double precision wmat(11,11),matvsk(10,10),matvku,covkusk(10),beta(3),m42&
&,m24,m222
integer i,j

forall(i=1:10,j=1:10)
	matvsk(i,j)=0.0d0
end forall
forall(i=1:11,j=1:11)
	wmat(i,j)=0.0d0
end forall

!!!
m42=3.0d0			!E(e1^4*ej^2)
m222=1.0d0			!E(e1^2*ej^2*ek^2)
m24=3.0d0    !E(e1^2*ej^4)
matvsk(1,1)=15.0d0
matvsk(6,1)=m42
matvsk(1,6)=m42
matvsk(8,1)=m42
matvsk(1,8)=m42

matvsk(2,2)=15.0d0
matvsk(4,2)=m24
matvsk(2,4)=m24
matvsk(9,2)=3.0d0
matvsk(2,9)=matvsk(9,2)

matvsk(3,3)=15.0d0
matvsk(5,3)=m24
matvsk(3,5)=m24
matvsk(7,3)=3.0d0
matvsk(3,7)=matvsk(7,3)

matvsk(4,4)=m42
matvsk(9,4)=m222
matvsk(4,9)=m222

matvsk(5,5)=m42
matvsk(7,5)=m222
matvsk(5,7)=m222

matvsk(6,6)=m24
matvsk(8,6)=m222
matvsk(6,8)=m222

matvsk(7,7)=3.0d0
matvsk(8,8)=m24
matvsk(9,9)=3.0d0
matvsk(10,10)=m222
!!!
wmat(1,1)=720
wmat(2:11,2:11)=matvsk
end subroutine var0




subroutine var0u(wmat,dmut,dvsig,pth,sigma,b)
!Variance of the moment conditions under the null
use ghs
use testmomsrot
use linear_operators
use matrixfuns
use datahadler
integer, parameter:: n=3,p=10
integer pth
double precision wmat(p+1,p+1),wmatmc(11,11),wmatu(pth+p+1,pth+p+1),dvsig(pth,n*n),dmut(pth,n)
double precision:: sigma(n,n),isigma(n,n),rsigma(n,n),irsigma(n,n),dmach,tol,b(n),beta(n)&
&,matV(n,n),irsigmah(n,n)&
&,matthmsk(n,p),matD(p+1,pth+p+1)
integer i,j,irank
!!!

isigma=.i.sigma
tol = 100.0*dmach(4)
call dchfac(n,sigma,n,tol,irank,rsigma,n)
rsigma=.t.rsigma
irsigma=.i.rsigma
beta=matmul(.t.rsigma,b)

if (dot_product(beta(2:n),beta(2:n))<1.0d-10) then
	matV=eye(n)
else
!	Sigma=eye(n)
   matV(1,:)=(/-0.57735026918963d0, 0.81649658092773d0,0.0d0/)
   matv(2,:)=(/-0.57735026918963d0,-0.40824829046386d0,-0.70710678118655d0/)
   matv(3,:)=(/-0.57735026918963d0,-0.40824829046386d0, 0.70710678118655d0/)

!	matV(:,1)=beta/dsqrt(dot_product(beta,beta))
!	matV(1,2:3)=(/-0.50497514502648d0,-0.58533678384716d0/)
!	matV(2,2:3)=(/0.78236638609871d0,-0.0d0/)
!	matV(3,2:3)=(/-0.36455855607615d0,0.81079026232156d0/)

!	print*, "Please, supply orthogonal matrix!"
!	print*, "beta", beta
!	stop
	!matV(:,1)=beta/dsqrt(dot_product(beta,beta))
end if

irsigmah=matV.tx.irsigma


call var0(wmatu(pth+1:pth+p+1,pth+1:pth+p+1))

forall(i=1:n,j=1:p)
	matthmsk(i,j)=0.0d0
end forall
matthmsk(1,1)=3.0d0
matthmsk(2,2)=3.0d0
matthmsk(3,3)=matthmsk(2,2)
matthmsk(2,4)=1.0d0
matthmsk(3,5)=matthmsk(2,4)
matthmsk(1,6)=matthmsk(2,4)
matthmsk(3,7)=1.0d0
matthmsk(1,8)=matthmsk(2,4)
matthmsk(2,9)=1.0d0

wmatu(1:pth,1:pth)=matmul(dmut,matmul(isigma,.t.dmut))&
&+.5d0*matmul(dvsig,matmul(kron(isigma,isigma),.t.dvsig))

wmatu(1:pth,pth+1)=2.0d0*(n+2.0d0)*matmul(dvsig,vec2(isigma))
wmatu(1:pth,pth+2:pth+p+1)=matmul(dmut.xt.irsigmah,matthmsk)
forall(i=1:pth,j=pth+1:pth+p+1)
wmatu(j,i)=wmatu(i,j)
end forall

matD(:,pth+1:pth+p+1)=eye(p+1)
matD(:,1:pth)=-matmul(wmatu(pth+1:pth+p+1,1:pth),.i.wmatu(1:pth,1:pth))

wmat=(matD.x.wmatu).xt.matD
end subroutine var0u


subroutine cov0(matcov,dmut,dvsig,pth,sigma,b)
use ghs
use testmomsrot
use linear_operators
use matrixfuns
use datahadler
integer, parameter:: n=3,p=10
integer pth
double precision matcov(pth,p+1),dvsig(pth,n*n),dmut(pth,n)
double precision:: sigma(n,n),isigma(n,n),rsigma(n,n),irsigma(n,n),dmach,tol&
&,matthmsk(n,p),b(n),beta(n),matV(n,n),irsigmah(n,n)
integer i,j,irank
!!!

isigma=.i.sigma
tol = 100.0*dmach(4)
call dchfac(n,sigma,n,tol,irank,rsigma,n)
rsigma=.t.rsigma
irsigma=.i.rsigma
beta=matmul(.t.rsigma,b)

if (dot_product(beta(2:n),beta(2:n))<1.0d-10) then
	matV=eye(n)
else
!	Sigma=eye(n)
   matV(1,:)=(/-0.57735026918963d0, 0.81649658092773d0,0.0d0/)
   matv(2,:)=(/-0.57735026918963d0,-0.40824829046386d0,-0.70710678118655d0/)
   matv(3,:)=(/-0.57735026918963d0,-0.40824829046386d0, 0.70710678118655d0/)

!	matV(:,1)=beta/dsqrt(dot_product(beta,beta))
!	matV(1,2:3)=(/-0.50497514502648d0,-0.58533678384716d0/)
!	matV(2,2:3)=(/0.78236638609871d0,-0.0d0/)
!	matV(3,2:3)=(/-0.36455855607615d0,0.81079026232156d0/)

!	print*, "Please, supply orthogonal matrix!"
!	print*, "beta", beta
!	stop
	!matV(:,1)=beta/dsqrt(dot_product(beta,beta))
end if
irsigmah=matV.tx.irsigma

forall(i=1:n,j=1:p)
	matthmsk(i,j)=0.0d0
end forall
matthmsk(1,1)=3.0d0
matthmsk(2,2)=3.0d0
matthmsk(3,3)=matthmsk(2,2)
matthmsk(2,4)=1.0d0
matthmsk(3,5)=matthmsk(2,4)
matthmsk(1,6)=matthmsk(2,4)
matthmsk(3,7)=1.0d0
matthmsk(1,8)=matthmsk(2,4)
matthmsk(2,9)=1.0d0



matcov(:,1)=2.0d0*(n+2.0d0)*matmul(dvsig,vec2(isigma))
matcov(:,2:p+1)=matmul(dmut.xt.irsigmah,matthmsk)
end subroutine cov0


end module moments3